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SUMMARY 


A fast and versatile procedure for algebraically generating boundaiy 
conforming computational grids for use with finite-volume Euler flow solvers is 
presented. A semi-analytic homotopic procedure is used to generate the grids. 
Grids generated in two-dimensional planes are stacked to produce quasi-three- 
dimensional grid systems. The body surface and outer boundary are described 
in terms of surface parameters. An interpolation scheme is used to blend 
between the body surface and the outer boundary in order to determine the field 
points. The method, albeit developed for analytically generated body 
geometries is equally applicable to other classes of geometries. The method 
can be used for both internal and external flow configurations, the only 
constraint being that the body geometries be specified in two-dimensional 
cross-sections stationed along the longitudinal axis of the configuration. 
Techniques for controlling various grid parameters, e.g., clustering and 
orthogonality are described. Techniques for treating problems arising in 
algebraic grid generation for geometries with sharp corners are addressed. A 
set of representative grid systems generated by this method is included. 
Results of flow computations using these grids are presented for validation of 
the effectiveness of the method. 

Introduction 

Numerical solution of field problems in any discipline, e.g., fluid dynamics, heat 
transfer, structures etc., requires a spatial discretization of the physical domain; 
that is, the generation of an organized set of points on the domain of interest for 
subsequent numerical solution of the governing equations. In the recent past 
numerous methods for generating grid-systems have been developed. These 
procedures compete with each other in their relative abilities to produce 
improved grid qualities such as smoothness of point distribution, clustering of 
points in regions of physical interest, ability to accommodate arbitrary boundary 
geometries, etc. A dominating philosophy in the realm of grid generation has 
been the introduction of a boundary conforming coordinate system[1] wherein 
coordinate lines are made coincident with physical boundaries thereby 
simplifying the application of finite - difference methods and eliminating the 
need for special procedures at the boundaries. Such systems allow 
transformation of the complex physical region into a uniform square mesh which 
simplifies coding at the cost of transformation of the governing equations 

The methods for grid generation in current use fall chiefly in two broad 
categories: those employing partial differential equations and those 
constructing grids by algebraic means. The method to be described here is an 
algebraic procedure to be used in generating grids for finite-volume methods. It 
should be made clear at the outset that the present grid system is boundary- 
fitted only in that after spatial discretization segments of the physical boundaries 
form sides of the computational cells on the boundaries. The finite-volume 
method diseretizes the computational region into hexahedrons bounded by 


quadrilateral planes. Thus, elemental volumes are bounded by flux surfaces 
and adjacent volumes share flux surfaces. The flux surfaces and volumes are 
computed from the actual physical coordinates of their vertices and therefore, 
there is no need for the introduction of boundary conforming curvilinear 
coordinate systems, the rigor of transformation into a uniform cube and the 
accompanying transformations of the governing equations. 

Algebraic grid generation procedures include such techniques as normalization 
of boundary curves, transfinite interpolation and the use of intermediate 
interpolating surfaces. An overview of algebraic techniques in current use has 
been given by Thompson in [2]. Algebraic techniques are attractive because 
they are fast, easy to implement and require no numerical solution of partial 
differential equations. On the other hand, algebraic methods generally lack the 
smoothness inherent in elliptic grid generations and are more sensitive to 
specified point distributions on the boundaries. Elliptic methods, however, are 
not entirely free from this sensitivity but are generally more tolerant of 
nonsmooth boundary data than algebraic methods. A method for preventing 
effects of boundary discontinuities from propagating infield will be described in 
a later section. 

Methods for grid generation by algebraic transformations fall in two broad 
categories: (1) shearing transformations and (2) boundary value 

interpolations. The present method is closely related to the second class of 
methods. In methods using shearing transformations, the four sided region 
enclosed by two curves connecting two straight lines in the physical plane is 
transformed into a rectangle by simply normalizing the coordinates. The points 
are distributed along the straight boundaries in the transformed rectangle 
according to a specified stretching function in order to obtain desired 
concentration of points. Clustering of points on the boundary is necessary in 
nonsmooth regions to prevent inward propagation of this lack of smoothness. In 
the boundary value interpolation technique a function on a region is found by 
interpolation between functions specified on the region boundaries. This 
technique is generally termed transfinite interpolation; blending function and 
connecting function interpolations being synonyms for the method. This 
procedure was used by Smith[3] to generate grids for wing-fuselage 
combinations using linear blending functions. A detailed description of 
transfinite interpolation has been given by Gordon and Hall[4]. Parameterization 
of the boundary curves should be carefully specified in order to avoid 
overlapping of grid lines. This method is more sensitive than partial differential 
equation methods to boundary point distributions. The blending functions may 
be polynomials whose arguments are expressed in terms of the curvilinear 
coordinates through stretching functions. An alternative way is to design the 
blending functions themselves to achieve the desired concentration of lines as 
in Eriksson[5] and Rizzi and Eriksson[6], where exponential functions are used. 
Grids generated algebraically in two-dimensional planes can be stacked to 
produce three-dimensional systems as in the present method. Such stackings 
may not in all cases produce sufficient smoothness in the direction of stacking. 
The present method provides for three-dimensional interpolation if desired. 
Application of transfinite interpolation in three-dimensions has been given by 
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Anderson and Spratley[7]. Eiseman[8] has applied the interpolation to a vector 
field of tangents in a multisurface transformation instead of the coordinate 
vectors. The tangents are defined on successive selected surfaces by the vector 
differences between positions. The interpolant creates a smooth vector field of 
tangents to the lines connecting the successive selected surfaces and the 
boundaries. The multisurface method has been extended in Eiseman and 
Smith[9] to use piecewise linear local interpolants in place of polynomials to 
provide local control over the coordinates. The controls may be used for 
example, to restrict the inward propagation of boundaiy slope discontinuities or 
to embed general rectilinear Cartesian systems within a global system while 
maintaining continuous first derivatives. 

The present method generates a quasi-three-dimensional grid system by 
means of a homotopic blending between inner and outer surfaces. The method 
was originally developed for blended wing-body geometries that were 
generated by a semi-analytic procedure. A method was developed by 
Barger[10] for designing blended wing-body configurations with cross-sections 
that vary smoothly from a initial prescribed nose shape to a prescribed base 
shape such that the cross-sectional areas conform to a preassigned axial area 
distribution. The two end cross-sections are specified by the designer and the 
intermediate cross-sections are computed at specified stations along the axis of 
the body by interpolation using blending functions. Grids are then generated in 
two-dimensional planes located at each of these stations. The method in its 
essence consists of finding a family of transition curves by blending between 
two parametrically represented boundaries. The inner boundary is the body 
surface and the outer boundary is specified according to the needs of the 
particular problem. If the configuration resides within a duct, e.g., a wind tunnel 
test section, the duct wall serves as the outer boundary, otherwise a circular or 
elliptic outer boundary is chosen. The grid points in space are obtained by 
interpolation and blending functions which depend on the boundary data. The 
scheme is devised from an extension of the homotopic procedure used to 
design the body surface itself. This basic procedure has been previously 
applied by Moitra[1 1 ,12] to generate grids for blended wing-body geometries. 
The two dimensional grids are finally connected in the body axis direction. The 
outer boundary shape is represented as a transition surface with the same 
distribution of parameter values as are used to describe the inner surface. In 
case of a nonanalytically generated body surface the discrete points can be 
parameterized with respect to some chosen quantity such as angular location, 
arc length, etc, thus making the method applicable to a general class of body 
geometries. 

The application of the homotopic blending method to grid-generation was 
treated quite briefly in reference 10. Only one example was given for a quasi- 
three-dimensional grid, and one for a fully three-dimensional grid. In view of 
this cursory treatment the purpose of the present paper is to pursue this subject 
in a greater depth and, in particular, to study the manner in which grid 
orthogonality at boundaries can be treated in the context of the homotopic 
blending theory. 
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Theory And Mathematical Development 

Algebraic schemes are known to provide the fastest and the least expensive 
ways to generate grid systems. The method used in this study is an algebraic 
scheme devised to generate body geometries as well as the associated grid 
systems which has been applied to blended wing-body geometries, fighter 
aircraft configurations and certain internal flow geometries. The method for grid 
generation is a natural extension of that for defining the body surface and will 
be better understood if viewed in conjunction with the body geometry 
generation scheme. 

Surface Equations: The problem of analytically designing a supersonic 
airplane or missile is generally influenced by two factors. One of these is the 
axial area distribution and the other is the cross-sectional shape, which except 
for axisymmetric missile bodies, varies from the nose to the base. Near the nose 
the cross-section is very often circular or nearly circular. Progressing along the 
longitudinal axis, the cross-section acquires various shapes as additional 
components such as canopies, wings, stabilizing surfaces are encountered. 
The same considerations as to shape variations are encountered in ducts and 
nozzles as well, since the inlet and exhaust cross-sectional shapes are usually 
different in realistic channel configurations. 


A method has been developed[10] for designing body geometries with cross- 
sections that vary smoothly from an initial prescribed nose shape to a 
prescribed base shape such that the cross-sectional areas conform to a pre- 
assigned axial area distribution. The essential elements of the method are 
presented in the present paper. The body geometry may be defined as a 
surface component such that it varies gradually and smoothly between 
specified initial and terminal cross-sectional shapes. The surface component 
may represent a forebody, a wing-fuselage combination or a duct. The direction 
of variation is chosen to be the x-direction, i.e., the longitudinal axis of the 
configuration and the two end shapes are defined in y-z planes. The surface is 
generated by devising a way to vary the size as well as shape of the 
intermediate sections lying between the specified end shapes. The initial and 
terminal cross-sectional shapes are assumed to be represented in parametric 
form in terms of a parameter t. Let (y B1 , z B1 ) and (y B2 , z g2 ) denote the parametric 

representations of the initial and terminal shapes respectively, so that 


at x = x 1 


and at x = x 2 


^bi ” VbiW 

Z B1 = Z Bl(*) 


«a * » * *b 


^B2 ” VbzW 

Z B2 = Z B2^) 


(la) 


(l.b) 
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( 7 ) 

Zbi = f ,(<p)sin<p 
0 < cp < 2n 

with similar expressions for y B2 and z B2 . A variation of the transition surface is 
given by 


/ v . 0 - C> C 
y(x,cp) = ty B1 .y B2 

/ v , (i - C) C 
z(x,q>) - ^ Z B1 - Z B2 


( 8 ) 


where C and X are functions of x, and y B1 , z B1 , y B2 , z B2 are functions of (p. The 
transition surface in cylindrical coordinate is given by 

x(x,cp) = X(x)[Ri(<p)J .[Rd<p)\ (9) 

assuming that R 1 and R 2 are the polar representation of the end sections. 

The technique gives the user a large degree of control over the particular 
geometry to be generated. A base shape with sharp wing-body junctions 
causes no difficulty in generating a smoothly varying intermediate geometry. A 

judiciously chosen distribution of the parameter q> in any cross-sectional plane 
allows circumferential clustering of points in regions of geometrical interest, 
such as sharp corners, thus ensuring proper definition of complex body 
geometries, while a properly chosen transition function C(x) serves the same 
purpose in the axial direction. An appropriate distribution of x values serves to 
cluster grid planes along the axis of the configuration. This procedure 
guarantees a smooth transition from the initial to the base shape while 
maintaining a specified area distribution. The two end shapes are specified in a 
fashion that ensures a natural one-to-one correspondence with respect to the t- 
parameter values. This correspondence is valid at every intermediate section 
throughout the transition and precludes the possibility of t = constant lines 
developing waviness or crossing each other. The surface definition 
equations(eqs (5) - (9)) were developed by Barger[10] for numerically 
generating aerodynamic surfaces. Once the end shapes are specified by the 
designer, the intermediate sections may be computed and displayed 

graphically. The shape and scale functions C(x) and X can be varied 
interactively until a desired transition surface is obtained. A schematic 
representation of the body definition procedure is presented in Figure 1. The 
prescribed shapes were a smooth circular shape at the nose SI, and a base 
shape S2 containing the central body and the trailing edges of the wing of a 
blended wing-body configuration. All the intermediate shapes, S3, S4, etc., 
were generated by blending between SI and S2. 
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It should be noted that if the end shapes are difficult to describe analytically it 
may be necessary to represent them numerically by means of digitized data. For 
those cases the data may need to be smoothed by b-splines. Similarly, the 

functions controlling shape and size variations C(x) and X respectively may be 
numerically specified; however, since these are normally smooth variations for 
aerodynamic shapes of interest, analytic expressions usually suffice for their 
representation. 


Grid Generation: A simple extension of the homotopic procedure for 

generating the body surface itself results in a scheme for gridding the external 
field. The (x,t) coordinate lines defining the body surface form a natural surface 
grid. For calculation of supersonic flow over such a surface a quasi-three- 
dimensional grid system is sufficient. The external field-grid is generated by 
finding a family of transitional curves between the body surface as the inner 
boundary and a prescribed outer boundary. The outer boundary is represented 
as a transition surface with the same distribution of x and t values as are used to 
describe the inner surface. In case of nonanalytic body geometries the discrete 
data defining the surface may be parameterized with respect to some 
convenient quantity such as angular location, arc length, etc. in each two- 
dimensional cross-section. The grid is generated by defining at each x-station a 
homotopic variation from the inner surface S^.z,) to the outer boundary 

S 0 (y 0 ,2 0 )- Denoting this homotopic parameter by r\, and the shape transition 
function between the boundaries by C*, we have, 


C to) 
C fo) 


0, T1 = 0 

1, n = i 


( 10 ) 


This function is taken to be independent of x; it retains the same value at each x- 
station. Essentially, the family of transition curves define a family of maps given 
by the homotopy 


{h^: A->B|r|el} (11) 

where I is the unit interval [0,1]. The inner and the outer boundaries Sj and S Q 
are two homotopic maps such that 

h o = 8, 

and (12) 

h , - s„ 

Next, the size variation function for the grid lines is denoted by X* so that 


X t < X < A, 0 


( 13 ) 
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where X. and X Q are the scaling functions used to generate the inner and the 
outer boundaries respectively. Since the size of the inner surface and possibly 
the outer surface vary with x, X is not independent of x. It is expressed in terms 

of the normalized function X by 


X = Xj + A.( X Q — A,j) 

so that 


X(T]) 

= o, 

r| = 

m 

= i, 

T1 = 


and X is independent of x. 


(14) 


(15) 


The grid is defined at each x-station in terms of the function C (-q) and M'h) ) and 
a natural correspondence is established between the grid points at the various 
x-stations. Smooth variations of the inner and outer boundaries with respect to x 
will result in smoothly varying grid lines in the x-direction. Appropriate choices 

for the function C (rj) and the distribution of rj result in approximate 
orthogonality of grid lines at the physical boundaries and clustering of grid lines 
near boundaries. The body-definition scheme is logically suggestive of a 
possible procedure for grid generation based on similar considerations of 
blending between an inner and an outer boundary. Finding a family of transition 
curves spanning the region between the boundaries consequently becomes the 
ideal procedure for generating the computational grid. A detailed description of 
the procedure follows. 

Coordinates of each boundary are expressed in terms of a parameter x For the 
inner boundary 


y i = yiCO 

Zj = Z|(l) 


(16) 


while identical expressions denote the outer boundary. In polar coordinates x 
could represent the angular coordinate. If z can be expressed as a single- 
valued function of y, then x = y. For shapes of greater complexity the choice of 
variables will vary. The procedure for grid generation consists mainly of 
determining a family of curves representing a gradual transition from the given 
inner boundary ty.Zj) to the outer boundary(y o ,z o ) in two-dimensional planes at 
each x-station. Assuming that the body surface coordinates have been 
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generated, a distribution of r| values is specified. The redistribution may be 

specified by means of polynomials, exponents, trigonometric functions etc., 

while ensuring that r\ = 0 on the body surface curve and rj = 1 on the outer 

★ 

boundary. The shape transition function C for the grid is then specified by 

* m 

C (ri) = 1 - ti (17) 

where m is a positive exponent that provides for additional control over line 

spacing near boundaries. Thus C = 1 on the inner boundary and C = 0 on 
the outer boundary with smoothly varying values distributed between 
boundaries. An arithmetic averaging between the boundaries yields the 
simplest family of grid lines given by 

y(r|,T) = C (Ti)yj(x) + [1 - C (Ti)]y 0 (x) 

(18) 

Z(T|,x) = C (Tl)Zj(x) + [1 - C (n)]Z 0 (x) 


where 0 < r| < 1 


For complex boundaries, e.g, those involving corners, a geometric averaging 
yields smoother transition curves given by 

* ★ 


yft.*) 

Z(T|,X) 


c (Tl) [1 - C (T|) 
Vito -y oW 


, C (Tl) [1 - C (Tl)] 


(19) 


The similarity with the scheme for body surface generation is clearly evident. 

Essentially C (ti) serves a purpose identical to that of C(x) in the body- 

* 

generation scheme. Here C provides for smooth blending between the shapes 
of the inner and the outer boundaries. Particular constant values of the rj- 
parameter generate specific curves of the family described by equations (18) or 
(19) and a pre-chosen distribution of r| determines the spacing of the resulting 
set of curves in each cross-sectional plane. 

A major advantage of the present procedure derives from the fact that it allows 
an integration of the procedures for body and grid generation provided the body 
cross-sectional shapes are definable by analytical means. For body geometries 
available as discrete sets of points only the grid-generation part is of interest. If 

the parameters t and x are chosen to be identical, e.g., the angular coordinate in 
both cases, then the body and the grid are definable by the same set of 
equations. It should be noted that the choice of the interpolation scheme, i.e. , 
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arithmetic averaging or geometric averaging for body definition and grid 
generation can be made independently of each other even though the extent of 
complexity of the body geometry usually has some influence on the choice of 
the blending technique for the grid. The final equations for generating y and z 
values for the entire field are found by combining equation (5) or (8) with 
equation (18) or (19). The four possible sets of grid-generation equations are 
given below. 

1) Arithmetic averaging in body definition and geometric averaging in grid 
generation( equations (5) and (19)): 


y(x,x,C,C ,tj) (x)[\(x){[1 - C(x)]y B1 (x) + C(x)y B2 (x)}] 

l(1 - C*(T1)) 


C (Tl) 


. VoM] 


( 20 ) 


z(x,x,C,C ,T|) = X (x)[X(x){[1 - C(x)]z B1 (x) 

C*fn1 (1 - C (rj)) 

+ C(x)z B2 (x)}] U [z 0 (x,x)] 

2) Arithmetic averaging in body definition and arithmetic averaging in grid 
generation ( equations (5) and (18)): 


y(x,t,C,C ,T|) = X (x)[[X(x){[1 - C(x)]y B1 (x) + C(x)y B2 (t)}]C (ti) 

+ (yo(x.t)lll -c'(r,)]] 

2 (x,t,C,C*,t|) = X (x)[[X(x){[1 - C(x)]z B1 (x) + C(x)z B 2 (x)}]C*(n) 

+ [z 0 (x,x)][i-c‘(n)ll 


( 21 ) 


3) Geometric averaging in body definition and geometric averaging in grid 
generation (equations (8) and (19)): 


* ‘y yr (1 - C(x)) 

y(x,x,C,C ,tj) = X (x)[X,(x)y Bi (x) .y B2 (x 


> CW ! 


c (Tl) 


[yo(x^)] 


(i - o (Tl)) 


z(x,x,C,C ,rj) = X (x)[X.(x)z B i (x) 


(1 • C(x)) C(x1 
.ZB2(^) 1 


C (ti) 


( 22 ) 
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4) Geometric averaging in body definition and arithmetic averaging in grid 
generation (equations (8) and (18)): 

* * (1 - C(x)) C(x) * 

y(x,x,C,C ,T|) = X (x)[{X(x)y B1 (x) .y B2 W }C ft) 

+ [y 0 (x,x)[1 - C (ti)]]] 

(23) 

* * r (1 -C(x)) C(x) * 

z(x,x,C,C ,ri) = X (x)[{X(x)z B1 (x) V .z B2 (x) '}C ft) 

+ [y 0 (x,x)[i - c*ft)]]\ 

The prescribed quantities are the nose and base shapes (y B1 .z B1 ) and (y B2 ,z B2 ), 
respectively, and the outer boundary (y o ,z Q ). All sets of curves including the 

body surface are then generated by a single set of equations. Therefore, at any 
given x-station, the body geometry and all curves between the inner and the 
outer boundaries are computed by a single procedure. This procedure is 
completed for each x-station. Spacing of the grid lines in all three coordinate 

directions can be controlled by choosing appropriate distributions of x,C,C ,x, 

and t). A schematic of the resulting quasi-three-dimensional grid system is 
presented in Figure 2. 

Some flow situations might demand a truly three-dimensional grid and a quasi- 
three-dimensional system may not be adequate; for example, a modification of 
the basic grid generation scheme allows the development of a fully three- 
dimensional grid as follows. If the influence of the body shape is required to be 
felt upstream of the nose, the origin in the axial direction can be defined at an 
X! + x 2 

average x-station — ;> along the axis of the body. Consequently the surface 

is now defined on the domain 

- <Pi,i ^ <Pi,i (24) 

while the t-domain remains unchanged 

-1 < t < 1 

The normalized outer boundary shape is defined to be an appropriate convex 

surface on the same (p , t domain with identical distribution of the <p and t 

variables as were used for the inner boundary. This normalized surface r(<p,t) is 

multiplied by a large scale factor X Q to determine the final outer boundary. A grid 

* 

homotopy parameter r| and a shape transition function C (q) are defined as 



done previously. The size transition function is denoted by A. (r|) such that 1 
< A < A Additionally a homotopic function |3((p) is defined such that |3 varies 

monotonically from 0 to 1 as |<p| varies from 0 to 1 . In the resulting scheme the 
grid-point locations are given by the following set of equations. 


x(<p,Ti) = A (ji)p(cp) 

y(cp,t,n) = A (ri){C (ri)yi(cp,t) + [1 -C (r|)]y 0 (r|,t)} 

( 25 ) 

Z(cp,t,Tl) = A Cn){C (Tl)Zj((p,t) + [1 -C (Tl)]Z 0 (Tl,t)} 

where (y^Zj) and (y Q ,z 0 ) denote the inner and outer boundary points 
respectively. 

An alternative way is to input different distributions for the x-values for the inner 
and the outer boundaries. Let us denote these distributions by x.(I) and x Q (I) 

where I denotes the I** 1 cross-sectional station. At each I-station a scale factor 
A is defined such that 

O 


A 


O 


Xo 

Xj 


(26) 


* 

Then a distribution for size transition A x (I) is defined for x-values of grid points 
* 

so that 1 < A x < A 0 . A blending equation for the x-coordinates of grid-points is 

then included in the set of grid generation equations. The x-values are 
generated by 

xO.rO = X x (I,Ti){C*(Ti)x j (I) + [1 - cWoffi) < 27 > 

Similar to the quasi-three-dimensional case, independent manipulation of the 
shape and size variation functions and redistribution provide considerable 
control over grid characteristics such as orthogonality, clustering etc. 

Generation of a grid around even simple wing-body configurations is by no 
means a simple task. The difficulties inherent in obtaining numerical solutions 
on a space discretized by a grid results in very demanding requirements for the 
grid generator that is expected to provide an optimum resolution with a 
minimum number of node points. It has been aptly pointed out by Eriksson[13] 
that the two chief components of computational fluid dynamics: grid generation 
and flow solution, should ideally be distinct operations and should be treated in 
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independent and modular ways. The method for grid generation should in no 
way be dictated by the particularities of a procedure for flow solution, and the 
flow solver should be able to accept as input an arbitrary set of points defining a 
grid and carry out solution of the governing equations thereon. In cases of 
practical interest, however, these two operations are intrinsically connected 
through such factors as data-structure, location and shape of boundaries and 
topology of the computational space. The present method of constructing three- 
dimensional grids for complex configurations is a homotopic blending 
procedure conceptually akin to transfinite interpolation and is attuned to the 
needs of finite-volume Euler flow solvers. This method generates grids by direct 
algebraic interpolation and not through algebraic transformation which is the 
norm with most transfinite methods in current use. A finite volume flow solver 
requires the physical coordinates of the grid points as input and transformation 
of the physical space into a cube in computational space is not imperative as 
would be in the case of finite difference algorithms. In the present procedure 
stretching transformations are replaced by stretching functions embedded in the 
blending functions of the interpolation methods to provide control of the point 
spacing. This homotopic scheme is transfinite, transformational logistic 
differences notwithstanding, in that it involves interpolation between functions 
specified along curves or surfaces rather than among point values and 
therefore, matches the functions at a nondenumerable set of points, although in 
practical situations these functions may be defined by discrete sets of values. 

The present technique, being an algebraic procedure , is fast and allows 
explicit control of grid point-distribution. The approach is particularly well-suited 
for use in conjunction with interactive graphics. Human intervention in the 
control process of grid generation will remain indispensable until the advent of 
techniques that are truly coupled with the equations of motion and grid control 
is entirely adaptive. With interactive graphics this intervention can be 
instantaneous since algebraic grid generation methods are explicit and require 
relatively few computations. Grids may be generated, viewed and revised until a 
satisfactory grid is obtained in a reasonably short period of time. 


Orthogonality 


Introduction 

Orthogonality of coordinate lines is a very desirable property of any grid system 
to be used in solving partial differential equations. Strict orthogonality is not an 
absolute necessity although excessive departures from it will result in 
deterioration of the numerical solution. Substantial deviations from 
orthogonality, however, are allowed by various flow solver algorithms. It may be 
prudent to point out that in the framework of general transformed computational 
spaces and finite difference procedures the use of strictly orthogonal systems 
allows the simplification of the transformed partial differential equations by 
dropping the terms in derivative transformations arising out of nonorthogonality, 
but the solutions can be substantially contaminated if any departure from true 
orthogonality should appear anywhere in the grid. In contradistinction to this the 
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retention of the full transformed equations and the use of a nearly orthogonal 
grid system in finite difference schemes has the advantage of introducing no 
numerical errors due to lack of orthogonality. The present grid generation 
scheme was developed particularly for use with finite-volume solutions, and 
since finite-volume solutions in physical space can be shown to be equivalent 
to finite-difference solutions in transformed space of the fluid dynamic equations 
in strong conservation form [14], it follows that near orthogonality in physical 
space is sufficient. Therefore, the terms orthogonality and near-orthogonality 
are considered synonymous in the present report. 

Methods currently in use for generating nearly orthogonal grids are generally 
based on the principle of adjusting an existing non orthogonal grid. The 
methods use various approaches such as grid generation by partial differential 
equations, conformal transformations and construction of orthogonal 
trajectories. As mentioned before, the resulting systems are not sufficiently 
orthogonal to allow dropping of the terms arising out of nonorthogonality in 
derivative transformations. 

Near orthogonality in P.D.E. grids is usually achieved by determination of 
orthogonal control functions. Elliptic generating systems produce coordinates 
that are more nearly orthogonal than those produced by most algebraic 
systems. The generating system of equations most often used is the Poisson 
system given by 


^xx + Syy = PfeTl) 
Tlxx + ^lyy = Qft’l) 


( 28 ) 


where P and Q are functions used to provide grid control. Thompson, et al[15] 
incorporated exponential control functions in P which were used by Steger and 
Sorenson[16] to consider iterative determination of constants involved in order 
to achieve orthogonality at the boundaries. A different approach was taken by 
Visbal[17] to generate nearly orthogonal coordinate systems for simply 
connected regions employing the general generating equations for orthogonal 
systems by computing the control functions P and Q from a preliminary sparse 
grid. 


Attainment of orthogonality at the boundaries while preserving conformality is 
generally not considered a viable proposition. Moretti[18] has advocated the 
use of a sequence of elementary transformations e.g., Karman-Trefftz, 
Joukowski, bilinear, etc., in succession to map an arbitrary curve to a near 
circle. Finally a shearing transformation may be used to map the near circle into 
a rectangle. The grid so obtained is expected to be nearly orthogonal. 
Sockol[19] has used two conformal mappings followed by a numerical grid 
generation to generate nearly orthogonal grids for cascades. More recently, 
lves[20] has applied conformal techniques to a three-dimensional configuration 
involving an inlet center body combination. In this procedure a conformal 
mapping is used in two coordinate directions and the distribution in the third 
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direction is produced by algebraic stretching.The resulting three-dimensional 
grid is orthogonal in two directions. 

Methods based on the construction of orthogonal trajectories are particularly 
well suited to orthogonalization of an originally nonorthogonal grid. These 
methods have an advantage over methods designed to produce orthogonal 
systems directly in that they allow easier control over line distribution. Since one 
family of lines from an originally nonorthogonal system is retained, the methods 
provide greater flexibility in generating a final orthogonal system with desired 
line distribution. 

A procedure for construction of grids by orthogonal trajectories has been 
applied by Graves[21] and Graves and Hamilton[22]. In this procedure, curves 
of a nonorthogonal system are first generated by any chosen distribution 
function along straight lines connecting points on the inner and outer boundary 
curves. Orthogonal trejectories are then constructed by the following steps: (1) 
Intersection of the normals at selected points on the inner boundary with the 
next curve outward are located: (2) The normal directions at the intersection 
points are determined and new intersection points using these new normal 
directions are located; and (3) final points on the outer curve are then chosen to 
be midway between the corresponding pairs of intersection points. This 
procedure is applied successively between each pair of curves of the original 
system. The points determined by the above procedure are then connected by 
straight line segments to generate the final grid. The resulting lines will not be 
strictly orthogonal to either the inner or the outer boundary but the departure 
from orthogonality will be small. In a related procedure Ghia et al.[23], 
determined points on successive curves of a nonorthogonal grid such that the 
orthogonality condition 

9ij - 0 i * j 

is satisfied. This procedure produces grids qualitatively similar to grids 
generated by the method of Graves. The initial nonorthogonal grid may be 
generated in a variety of ways, however, algebraic transformations are 
preferred for this purpose over other means because of the explicit and 
continuous definition of the transformations provided by algebraic methods. 
Shearing, Hermite or multisurface transformations are adequate for this 
purpose. A natural correspondence between the orthogonal trajectories and 
images of characteristic curves under the given nonothogonal transformation 
has been given by Eiseman and Erlebacher[24], A general development of this 
topic has been presented by Eiseman[25]. 

It is generally held that algebraic procedures for grid generation, unlike P.D.E. 
systems have no inherent mechanism for producing orthogonal grids and 
simultaneously maintaining smoothness and even when methods, such as the 
ones described above, are used to achieve near-orthogonality, intersection of 
grid lines of the same family is an often unavoidable problem. In this paper a 
theoretical development and a method will be presented that exploits the 
pseudo differential nature of homotopic procedures to produce grids that mimic 
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P.D.E. grids with reasonable closeness so far as grid qualities e.g., 
orthogonality, smoothness and non-intersection of grid lines are concerned. 

Present Method: 

Understanding of the present method for achieving near orthogonality and 
smoothness from homotopic procedures will be facilitated by a review of the 
basic tenets of homotopy theory that play important roles in the development of 

k r + 1 

the theory. Given X c R Jet Xx[0,1] denote the subset of R consisting of 

all (x,t) with x e X and 0 < t < 1 . Then two mappings 


f, g : X -4 Y 

are called smoothly homotopic if there exists a smooth map 


(29) 


F : X x[0, 1] -» Y 

with 

F(x,o) = f(x) 

F(x,1) = g(x) 

where t denotes the homotopy parameter. 


(30) 


(31) 


Many of the properties of a map are not altered if the map is deformed in a 
smooth manner. Intuitively, a smooth map f : X -4 Y is a deformation of another 

f 0 : X -4 Y if they may be joined by a smoothly evolving family of maps.f t : X -» Y 
In the present method of grid construction equations (18) or (19) may be seen to 
be algebraic representations of the homotopic family f^ that smoothly maps the 

inner boundary into the outer boundary through the homotopy parameter ri that 
varies between 0 and 1 from one boundary to the other. Homotopy is an 
equivalence relation on smooth maps from x to y. In order to recover the 
smoothly evolving family of maps joining f and f 1 one defines. 


f t : X -4 Y by f t (x) = F(x,t) (32) 

where F(x,t) can be specified algebraically. 

It is assumed that f t is a diffeomorphism. A function f is a diffeomorphism if f 

exists and both Df and Df ' 1 exist and are continuous. This assumption is not 
restrictive and has the following consequences. 
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(i) f t (x) = f t (y) if and only if x = y, hence trajectories are uniquely 

specified by their initial conditions. 

(ii) The derivative of a trajectory with respect to the initial condition 

exists and is non singular. 

Consequence i) may be seen as a statement of one-to-one correspondence. 
Both consequences are automatically satisfied by choosing F to be an explicitly 
specified algebraic relation. 

It can be shown[26] that diffeomorphisms constitute a stable class of smooth 
maps of a compact manifold x into a manifold y .The present method of attaining 
orthogonality through homotopies is based on the concept of stability of maps. 
From a point of view of practicality and physical measurements no continuous 
functional relationship is ever perfectly determined over a wide range of 
argument values. Consequently, the properties of a mapping that are of 
physical relevance are only those which remain valid when the map is 
subjected to slight deformations. Such properties are termed "stable" properties, 
and the collection of maps that possess a particular stable property may be 
referred to as a stable class of maps. In essence, a property is stable provided 
that whenever f : x - y possesses the property and f : x - y is a homotopy of f , 

then for some e > 0, each f t with t < e also possesses the property. Thus small 
perturbations in the homotopy parameter do not alter the stable properties of a 
map. 

Here we are concerned with the orthogonality of the trajectories joining 
corresponding points on the successive curves, specially at the given 
boundaries. In this context a few stability considerations are particularly 
important. The properties of grid smoothness and conformity of the overall grid 
with the given boundaries are stable under slight deformations of the map. The 
property that a curve pass through a given point is not stable, since a small 
change in the homotopy parameter can immediately distort any such curve to 
avoid the point, (Fig. 3). It may be seen therefore, that the points on a curve may 
be redistributed by judicious choice of the homotopic parameter distribution in 
order to achieve a desired set of trajectories. The points on the boundaries are 
of course fixed and not altered by these changes. Another property of interest is 
that of intersection of curves. The naive point-set condition of intersection is 
seldom stable, and therefore, is meaningless in the physical world. 
Nontransversal intersection, however, is physically relevant in the present 
context , and this property can be readily seen (Fig. 4) to be unstable. Two 
curves of the family that intersect at an even number of points may be separated 
by adjusting the distribution of the homotopy parameter. These properties 
considered together suggest the possibility of generating approximate 
orthogonal trajectories while preventing intersection of curves of the same 
family by choosing appropriate distributions of the homotopy parameter. These 
particular distributions can be determined from the boundary conditions. In 
figures 5.a and 5.b a family of curves before and after small perturbations in the 
homotopy parameters for achieving near- rthogonality of the trajectories, are 
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presented schematically. The perturbations redistribute the points lying on each 
curve so that the trajectories joining them become orthogonal. 

The problem of intersecting trajectories is presented graphically in figure 6. For 
extreme cases of boundary curves with sufficient concavity, an attempt to 
orthogonalize the trajectories may result in intersecting trajectories. In Figure 
6.a, two such trajectories emanating from the concave region of the boundary 
are seen to intersect at a pair of points. In the enlarged view presented in Figure 
6.b the two trajectories are seen to intersect at the points A and B. The curve f b 

lying between the points of intersection is seen to locally turn back on itself. An 
arbitrary curve of the family, f , may now be imagined, which intersects the curve 

f^ at two points. The problem of intersection of the trajectories at points A and B 

may thus be translated to the intersection of two curves of the same family f at 
points P and Q as seen in Fig 6.c. As mentioned before, this type of intersection 
is an unstable property of the map. The homotopy parameter b may therefore be 
locally adjusted to force f b completely clear of f p thus removing the problem of 

intersecting trajectories. The present procedure for generating near-orthogonal 
grids while preventing grid line intersections does not necessitate generating a 
non-orthogonal grid first and then adjusting this preliminary grid. The controlling 
values and distributions of the homotopy parameter are derived from the 
boundary data and used directly to generate the final grid. This capability of 
generating orthogonal grids is inherent in the homotopy procedure in much the 
same way as the method of orthogonal control functions is inherent to the P.D.E 
grid generation systems. It may be noted that the relation of smooth homotopy is 
an equivalence relation and consequently transitive [27]. This method of 
constructing orthogonal grids by homotopy is therefore equally applicable to the 
multisurface procedure. The algebraic description of the homotopy method of 
orthogonal grid construction follows. 


The basic interpolation scheme, equations (18) may be re written as 


where 


x = Xj E + x o (1 - E) 

y = y s e + y o 0 - e) 


(33) 

(34) 


The subscripts i and o denote the inner and the outer boundaries respectively. 
Here E varies between 1 at the inner boundary and 0 at the outerboundary and 
plays a role analogous to that of a homotopy parameter. Modifications of E, 
therefore, causes slight deformations of a given map and may be used to 
achieve orthogonality. It can be seen that E serves as a parameter that 
propagates the influence of the inner boundary on the coordinates of any point 
in the field. In order to achieve orthogonality of the trajectories at the inner 
boundary the influence of the inner boundary on points near it needs to be 
increased. Here, this control is provided through the use of an exponent P. The 
coordinates x and y of any point in the field may then be modified by 
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i 



X 


( 35 ) 


= Xj E P + x Q (1 - E P ) 

y = y. E p + y 0 (1 - E p ) 

P <1.0 

The exponent P is not a constant and must be determined from the boundary 
data subject to the orthogonality condition. Two kinds of segments of the inner 
boundary that need to be considered separately are a) those that are nearly 
horizontal and b) those that are nearly vertical. For nearly horizontal segments, 
orthogonality at the boundary may be obtained by modifying the x-coordinate of 
the point on the curve lying next to the boundary as shown in Figure 7.a. In 
nearly vertical segments a modification of the y coordinate achieves the same 
effect. (Fig 7.b). The exponent P for any trajectory is determined from the 
condition that two vectors be orthogonal. One of these vectors is found by 
connecting the point(Xj,yj) on the inner boundary and the point (x,y) lying just off 

the boundary on the trajectory in question (Fig. 8). The other vector passes 
through the print (Xj,yj) and a point ( x,y ) on the line passing through (Xj.yj) and 
parallel to the line joining (x j+1 ,y j+1 ) and (x j _ 1 ,yj_-| ) on the inner boundary. 

These two vectors are denoted by A and B respectively. The method for 
evaluating P for the two kinds of inner boundary segments are presented below. 

a) Nearly horizontal segments 

As mentioned before only the x coordinate needs to be adjusted for trajectories 
emanating from inner boundary segments that are nearly horizontal. The 
orthogonality condition is that the dot product of the vectors A and B be zero, 
i.e., 

A.B = 0 (36) 

which translates into 


(x - ^ )(x - ) + (y - yj )(y - y 4 ) = 0 


(37) 


Substituting 


x = X, EP + (1 - E P )X 0 


(38) 


and y = yj E + (1 - E)y o 

one obtains 


(Xq - x,)(i - E p )( x - X, ) + (y 0 - Vj )(1 - E)(y - Vj ) (39) 
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Rearranging terms results in 


E p = 1 + A(1 - E) 

where 


( 40 ) 


(y 0 - yjXy - yd 

/ 

(x o - Xi)(x - Xi) 


(41) 


Finally solving for P one yields 

_ ln(1 + A( 1 -E)) 
In E 


(42) 


where E has the value corresponding to the homotopic curve lying next to the 
inner boundary. The value of P given by equation(42); will result in near 
orthogonality between the trajectory emanating from (Xj.yj) and the inner 

boundary. 

b) Nearly Vertical Segments 

For boundary segments that are nearly vertical it is necessary to adjust the y 
coordinate values for points on the trajectory. The values of x and y are then 
given by 


x = XjE + (1 - E)x q 
y = Yi E p + (1 - E p )y 0 

substitution of (43) in the orthogonality requirement 
A.B = 0 


(43) 


results in 


(E - 1)(Xj - x 0 )(x' - Xj) + (E p - 1)(y, - y 0 )(y' - y s ) 


(44) 


Rearrangement of terms yields 
E P = 1 + B(1 - E) 


(45) 
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where 


B = (Xj - Xq)(x -Xj) 

(yi-yo)(y -vi) ( 46 ) 

The final solution for P is 

p = ln(1 + (1 - E)B) 

In E (47) 

The value of P given by (47) will result in near orthogonality of trajectories 
emanating from nearly vertical segments of the inner boundary. A logical switch 
from scheme a to scheme b depending on the nature of the boundary segment 
can be easily implemented. 

The problem of grid lines or trajectories intersecting each other in the interior of 
the field is considered next. As seen above the use of the exponent P derived 
from the orthogonality condition ensures orthogonality at the physical boundary. 
Intersecting trajectories can be separated by modifying P. This may be achieved 
by the use of an exponent q such that 

x = XjE* 3 + x 0 (1-E^) q<1.0 (48) 

Raising P to the power q increases the value of the exponent of E, which in turn 
makes the grid lines follow the inner boundary less closely and reduces the 
chance of trajectory intersection. Using a constant value for q, however, will 
reduce orthogonality. In order to maintain orthogonality near the boundary the 
logical solution is to cause q to decay as one proceeds outwards from the inner 
boundary along trajectories. This is achieved by making q a function of E such 


that 

q = F(E) 


(49) 

where 

q = 1 

at the inner boundary 


and 

q = o 

at the outer boundary. 



The use of equation (48) with q defined as in (49) will preserve grid 
orthogonality in the region close to the boundary while preventing grid line 
intersection away from the boundary. 

To assist in following the train of logic involved in this procedure of successive 
exponentiation a formal logic proof of the concept follows. The region of interest 
is understood to be the region near the inner boundary where orthogonality is 
required. The predicates O and G represent grid orthogonality and grid line 
intersection respectively and P,Q and Q e denote the acts of using the exponent 
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P, raising P to a constant exponent q, and making q a function of E in that order. 
In the following table the procedures of exponentiation, their formal logic 
representation and comments explaining the consequences are given. 


Procedure 

Formal Logic 

Comments 

EP 

P=>0 (1) 

Raising E to the power of P as in 
equation (35) produces ortho- 
gonality near the boundary 
while trajectories may or may 
not intersect 

E P V 

/q * const 

Q =>'(0 a C) 

=> '0 v 'C 
(Demorgan’s Law) 
=> P v ( C (2) 

(from (1 )) 

Raising P to the constant power 

q will remove grid line intersec- 
tion at the cost of orthogonality, 

i.e., will prevent the construction 
of grids that are orthogonal and 
contain grid intersection. 

E P V 

/q-f(E) 

Q e =>PaQ 

=> P A ('P v 'C) 
from (2) 

=»(Pa P) v (P a 'C) 

(Laws of distributivity 
of conjunction with 
respect to disjunction) 

[(Pa'P) is an absurdity] 
=>Pa'C 

=>0a' c 

Use of q = f(E) preserves P near 

the boundary while maintaining 
the effects of Q in the inner field. 

Result (3) shows that the use of 
q = f(E) produces orthogonality 
and eliminates grid line cross- 
ings. 


The procedure described above produces smooth grids that are nearly 
orthogonal near the inner boundary and do not contain intersecting grid lines. 
The procedure is simple since all blending and control functions are specified 
explicitly. As opposed to other methods of grid construction by orthogonal 
trajectories, it is not necessary to generate a preliminary nonorthogonal grid 
and later adjust it. All exponents and coefficients are derived from the boundary 
data and applied explicitly to generate the final grid in one execution of the 
procedure. 
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Description of the algorithm 


The algorithm is structured to produce nearly orthogonal, smooth grids for any 
arbitrary three dimensional configuration. The configuration geometry is 
assumed to be available in the form of two-dimensional curves in planes 
stationed along the longitudinal axis of the body(see Fig. 1). Required inputs 
are the the shape functions defining these cross-sectional shapes, factors 
governing their size variations, the shape of the outer far-field boundary and an 
optional file containing the discrete point-set representation of the body 
geometry. Brief discussions on important parts of the algorithm follow. 

Body surface definition 

a) Analytic Case: 

Blending between two prescribed end shapes is performed to generate all 
intermediate sectional shapes. Equation (5) defines the transition surface 
shown schematically in Figure 9, in terms of the end shape functions and two 
coefficients, each defined as a function of a single independent variable. The 
two end shapes (YB1 ,ZB1 ) and (YB2.ZB2) depend only on t, whereas the shape 
variation function c and the size variation function X depend only on x. The 
normalized end shape cross-sections are specified by defining the functions 
Fj(t) and Gj(t), where 

F-j(t) = ZB1 (t) 

G-j (t) = YB1 (t) 

(50) 

F 2 (t) = ZB2(t) 

G 2 (t) = YB2(t) 

and t is a circumferential parameter with a chosen distribution, which varies 
from -1 to +1 from one wing tip to another. The functions Gj(t) are usually set to 

have the corresponding t-values, i.e. 

G lW = 1 (51) 

G 2 (t) = t 

while F^(t) and F 2 (t) can be specified by any chosen algebraic formulae. F 1 and 
F 2 are usually given separately for the upper and lower surfaces of the 

configuration; however, care must be exercised to ensure that they match at 
appropriate points such as wing tips and wing-body junctions. The normalized 
cross-section at the nose of the configuration is usually an ellipse or a circle as 
schematically shown in Figure lO.a. The terminal cross-section typically 
contains a circular central body and the trailing edges of the wings as 
represented in Figure lO.b. Once the end shapes have been determined, a 
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normalized sectional shape at station x along the longitudinal axis is computed 
by combining the F 1 and F2 shapes as follows, 


FI = C(x)F 1 + [1 - C(x)]F 2 (52) 

GI = C(x)G 1 + [1 - C(x)]G 2 

where C(x) is any function of x that equals 1 at x n and 0 at Xf. The distribution of 

C(x) determines the nature of the blending between the two terminal shapes. 
The scaled cross-section at station x is then determined by multiplying the 
coordinates (FI,GI) at x by a scale function SC(x), which is normally zero at the 
nose(x n ), and given a specified size at the base at Xf. The final coordinates of 

the configuration surface (YI,ZI) are then given by the following equations: 

YI = SC(x){C(x)G 1 + [ 1 - C(x)]G 2 } 

ZI = SC(x){C(x)F 1 + [ 1 - C(x)]F 2 } < 53 ) 

b) Discrete Case: 


In the event that the surface geometry is not defined by analytic means, the 
geometry may be input to the code as an ordered set of discrete points. The 
input file must contain the grid points defining the cross-section in each axial 
plane in succession. Moreover, the number of discrete points in each plane 
must remain constant in order to ensure continuous surface grid lines. 


A normalized shape for the outer boundary is first prescribed through the shape 
functions Fg(t) and G 3 (t). The function F3 is specified as an algebraic relation 

defining a desired curve, e.g., a circle, and G 3 (t) is set equal to t as in the case 

of the configurations end sections. The shape of the outer boundary may be 
varied from one x-station to another by making F 3 a function of x. A schematic 

representation of the outer boundary shape appears in Figure 11. Once the 
outer boundary has been specified at each station, a set of grid lines is 
calculated in each plane by blending between the inner and outer boundaries 
through the parameter ri using equations(18). A distribution of rj -values along 
trajectories is first chosen and then the distribution of the grid lines is 
determined by the function E(r|) where 

E(ti) = C*(r|) as defined in equation (17) , where 

E(t|) must be 1 when r\ = 0 and 0 for r\ = 1. In order to ensure near orthogonality 

at the boundary and nonintersection of trajectories, the coefficient P of E(ti) has 
to be determined from equations (42) or (47) depending on whether the inner 
boundary is locally nearly horizontal or vertical. For this purpose a test on the 
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local slope of the inner boundary is performed. If the absolute value of the slope 
exceeds a preset limit the surface is considered nearly vertical and nearly 
horizontal otherwise. The actual scaled size of the outer boundary is defined by 
the multiplication factor SCL3. Consequently, the set of transition grid lines from 
the inner surface(SC.FI.SC.GI) to the outer grid boundary. (SCL3.F 3 , 
SCL3.G3) is given by 

Z = E.SC.FI + (1 - E).SCL3.F3 

Y = E.SC.GI + (1 - E).SCL3.G3 (54) 

where Z,Y are the final coordinates of the points in the field. Thus the 
parameters controlling the grid are x,t and E whose distribution may be varied 
as desired. The shape of the body can be changed by varying the end-shapes 

F 1 (t)> G-j (t) 

F 2 (t), G 2 (t) 

and the shape transition and scale functions 

C(x) 

SC(x) 

while F 3 (t), G 3 (t) determine the shape of the grid. 

The generated grid points are saved in a file for graphics post-processing or for 
use by a flow-solver code. A detailed description of the computer code for grid 
generation can be found in a companion report[28]. 

Examples of generated body geometries and arid systems 

A representative spectrum of body geometries and coordinate grid systems 
produced by the theory and numerical procedures described in the previous 
sections will next be presented. The versatility of the method is manifested in 
the presented examples. The effects on the body surface shape and grid 
control, of the various parameters, coefficients and functions used in the theory 
are illustrated. Computed grid examples include those for analytically defined 
body geometries as well as ones for body geometries represented by discrete 
point sets. 

Body geometries: 

A set of four body geometries generated by the method of blending between 
prescribed nose and base shapes are presented. For each case the shape 
transition, the scaling factors and their effects on the resulting body shape are 
explained. All of these geometries were generated using equations(5) of the 
body definition scheme. 
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Body 1. A three-view representation of this configuration is given in Figure 
12. The frames are appropriately labeled with the coordinate directions. The xy 
frame contains the planform shape, the xz frame represents the side elevation 
and the front projection appears in the yz plane. The configuration was 
generated by specifying elliptic shapes for both the nose and base cross- 
sections. The final configuration has elliptic cross-sections at each x-station with 
the major and minor axes in a ratio of 1 :0.18. The shape transition function C(x) 
and the scaling function X(x) used for this configuration are defined as 


C(x) = 1.0 


MX) 


x(12.0 - x) 
20.0 


(55) 


Setting C(x) equal to the constant value of unity throughout has the effect of not 
altering the shapes of the intermediate cross-sections at all, so that if the end 
shapes are elliptic the intermediate cross-sections are strictly elliptic also. The 
specified function of the scaling factor produces the curved outer edges 
appearing in the planform and the side elevation. 

Body 2. This configuration represents a spiked fore body with small canard 
surfaces. The nose shape is specified to be a very small circle and the base 
shape contains a circular central body and the straight line edges of the canard 
surfaces. The intermediate cross-sections gradually change in shape from the 
nose to the base as determined by the shape transition function C(x) and the 
scale factor X(x). The function C(x) is given by 

C(X) = ( 56 ) 

where x n and x f are the x values at the nose and the base, respectively. 


Two different distributions are used for the scaling factor X(x), one for x values 
lying between X n and one-third of the length of the configuration and another for 
the rest of the length. The distributions are given by 


X(x) = M + (l m - 1)7^- x< x 1/3 

x 1/3 


(57) 


Mx) - ^1/3 + (^f ^1/3)^ 


X - X 


1/3 


Xf - X 


X > X 


1/3 


1/3 
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where x 1/3 is the x-value for the point lying at one-third the length from the nose, 

and Xjf X f and X 1/3 are the prescribed scaling factor values at the nose, the 
base and at x 1/3 respectively. 


Body 2. A blended wing-body configuration is shown in Figure 14. This 
geometry is representative of modern high-speed aerospace vehicles or 
missiles which are characterized by smooth blending of the wing into the 
fuselage. As in the previous case, the prescribed end shapes are a small circle 
at the nose and a circular fuselage and wing trailing edges at the base. In this 

case the shape transition function C(x) and the scaling factor A.(x) are given by 


C(x) 




1.7 


. . . __0.8 . ^ —0.8. . 
X(x) = R A, + (1 - R A n 


(58) 


where H ~ Y . Y 

Af a n 

The value of the function C(x) varies from 1 at the nose to 0 at the base and 
generates the intermediate cross-section shapes by a smooth blending 
between the two end shapes. 


Body 4. This case illustrates the use of the scaling function X(x) to generate 
body shapes of varying complexity. The configuration presented in Figure 15 
has the same nose and base shapes as the previous case but the scaling 
function has been modified to produce the interesting planform of the blended 
wing-body geometry. The shape transition function C(x) is identical as that used 
for body 3. The shape controlling functions for this configuration are 


and 


where 


C(x) 

Mx) 

Q = 


= (^) 
x f - x„' 


1.7 


= Q + 0.15 Sin 


<^)K 


- A* n 

( 7~1T )(X ‘ x "> 

Xf - x n 


(59) 

(60) 


(61) 


and x 1/3 has the same definition as given for body 2. The intermediate cross- 
sections are produced by a smooth blending between the end shapes while 
conforming to the size variation imposed by equation (60). 
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These examples demonstrate that independent control over the functions 
controlling the shape and size variations can be used to generate complex 
aerospace vehicle shapes for use in high-speed computations. The technique 
is particularly suitable for generating advanced body shapes having smooth 
wing-body junctions, i.e. blended wing-body geometries, that eliminate the 
strong shocks that would occur at sharp wing-body junctures. 

Computed arid examples. 

The homotopic scheme is applied to several examples in the present section. 
Grid configurations obtained using the basic scheme are presented first. The 
techniques for orthogonality and prevention of grid intersection are illustrated in 
the examples that follow. The body geometries for which grids were generated 
vary in complexity ranging from simple airfoils to blended wing-body geometries 
and fighter aircraft. For most cases only a part of the field near the inner body 
surface is presented for clarity. Relevant values of coefficients and parameter 
distributions are given for each example. 

Basic grids. 

The grid examples in this group, presented in Figures 16 to 25 were generated 
using equations 18. No special provisions for orthogonalization or prevention of 
crossing of grid lines were incorporated in these cases. The lines of the 
homotopy map were uniformly distributed between the inner and outer 
boundaries through the following distribution functions 

K - 1 

11 N - 1 (62) 

and C (rj) = 1 - r| (63) 

where ri is the homotopy parameter, N denotes the index of the outer boundary 
and K varies between 1 and N from the inner to the outer boundary. 

The sectional planar grid shown in Figure 16, was generated for the elliptic 
body presented in figure 12. The outer boundary was prescribed to be flatter at 
the bottom in order to facilitate computation of high angle-of-attack flows. The 
smoothness of the transition of the lines of the map between the boundaries is 
apparent. The trajectories in this case are straight lines joining corresponding 
points on the inner and outer boundaries. Orthogonality was not imposed; 
however, the grid was found to be satisfactory for finite-volume computations at 
moderate Mach numbers. The orthogonalized grid for this configuration will be 
presented in a later set of results in this section. Sectional grids at two x-stations 
of the body geometry presented in Figure 14 are shown in figures 17 a and 17 
b. The sections presented in figures 17 a and 17 b are located at the midlength 
of the body and near the base, respectively. A smooth family of curves is seen to 
have evolved in each case; moreover the discontinuity at the sharp wing tips 
does not propagate into the field as would be usually expected for algebraic 
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procedures. Similar plots of sectional grids for two sections of the body shown 
in Figure 15 are presented in Figures 18-a and 18-b. 

The grid examples presented so far are illustrations of O-grid topologies in 
planes normal to the axis of the body and the direction of the flow. The next two 
examples are for C-grid topologies where the plane of the grid is parallel to the 
direction of the flow. A planer grid for a section of a semi-elliptic cylinder, 
typically representing the leading edge of a wing, is presented in Figure 19. The 
outer boundary is prescribed as a semi-circle of a chosen radius and the 
distribution of the family of curves is again specified to be uniform. A grid in a 
chord-wise plane of a supersonic wing is presented in Figure 20. The inner 
boundary in this case was specified by a discrete point set. The figure 
represents the region close to the airfoil and a smooth blending of the grid lines 
is evident. The grid lines have been distributed uniformly. 

The quasi-three-dimensional grid generation technique is illustrated in Figures 
21 a through 21 e for a typical fighter aircraft configuration. A planer grid for a 
section of the configuration containing the fuselage and a vertical stabilizer 
surface appears in Figure 21. a. Quasi-three-dimensional grid systems are 
constructed by generating a planer grid for each cross-section of the body and 
finally connecting the corresponding grid points on these planes in the body- 
axis direction producing a grid-system that completely envelopes the body. A 
cut-away view of the total grid system is shown in Figure 21 .b. A portion of the 
outer boundary surface has been cut away to reveal the underlying grid 
structure. The grid system is presented for the left half of the configuration and 
five planer grids at different locations along the axis of the configuration are 
visible. The configuration with two planer grids appears in Figure 21 .c. Enlarged 
views of these planer sections are presented in Figures 21. d and 21. e. 

The next set of results will demonstrate the effects of grid-stretching, i.e., 
nonuniform distribution of the grid lines through control of the homotopic 

parameter r\, and grid-clustering on the body surface by chosen distributions of 
the circumferential parameter t. Concentration of grid lines near the inner 
boundary is effected by prescribing values for the coefficient m in the definition 

of C*(r|) in equation (17). It has been shown before that a value of unity 
generates uniformly distributed lines. The planar grids presented in figures 22. a 
and 22. b show grid line concentrations achieved by setting m = 1.8 and m = 
2.2, respectively. Greater concentration of lines in the wall boundary regions is 
achieved by larger values of m. Similar plots for two sections of a blended wing- 
body configuration appear in Figures 23.a and 23. b. The chosen value for m is 
2.0 for both sections, however, m may be specified as a function of x in order to 
achieve varying degrees of concentration in different sections. A comparison of 
Figure 23 with Figure 18 reveals the additional grid concentration near the inner 
surface. A sectional grid for an airfoil wing section appearing in Figure 24 
shows grid lines concentrated near the wing surface using m=1.8. 

Grid clustering at any location on the body surface can be achieved by 
prescribing a modified cosine distribution for the circumferential parameter t. In 
a given section t varies from -1 to +1 from one wing-tip to another. For bodies 
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with a plane of symmetry it is sufficient to specify the t-distribution on only one 
side of the plane of symmetry from this plane to the wing-tip, i.e., in the domain 

0 < t < 1 The following examples illustrate grid-clustering at various points on 

the body surface along with the t-distributions that generated the clustering. A 
sectional grid with grid lines clustered at the tips of a body with elliptic cross- 
section is presented in Figure 25. The distribution function for t in this case is 
given by, 


A 


-cos 


(J - 1)ti 

TN^ry 


(63) 


t 


a (1 ' A) 


(64) 


where N is the total number of circumferential points on the surface in a given 
section and J varies from 1 to N circumferentially along the surface from one 
wing-tip to another. A similar distribution was used to cluster grid points at the 
wingtips of a blended wing-body configuration as shown in Figure 26. A 
modification of the exponent in equation (64) results in clustering near the plane 
of symmetry. The grid presented in Figure 27, showing clustering at the plane of 
symmetry was generated using 

, = A (i.o + F(1.0 - A)) (65) 


where A is given by equation (63) and the value of F determines the degree of 
clustering. Grid lines can be clustered near any point on the body surface 
denoted by t = t c by modifying the point distribution by the following equations. 


and 


[ P + Q ( 1-10 
t = A 1 l 1.0-tc 


)] 


t = A 


[R+S(!fi)l 


( 66 ) 


(67) 


where the constants P, Q, R and S are carefully chosen to provide desired 
clustering. The grid shown in Figure 28 was obtained by clustering neart c = 0.4 

for a section of a blended wing-body configuration. The values used for P,Q,R 
and S were 1.1, 1.2, 0.98 and 0.25 respectively. 

The grid examples presented so far were produced with no efforts at 
orthogonalization. The following group of examples demonstrate the imposition 
of near orthogonality at the inner boundary by the use of equations (38) and 
(43) with the values of the exponent P given by equations (42) and (47). A 
sectional grid for a body with elliptic cross-section, with partial orthogonality of 
grid lines near the body surface is presented in Figure 29. An enlarged view of 
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the region near the inner surface appears in Figure 30. Grid lines are seen to 
start nearly orthogonal to the body surface and then gradually curve away 
toward their destination points on the outer boundary. Orthogonality of the 
trajectories is more clearly visible in the sectional grid for a blended-wing-body 
configuration presented in Figure 31. The scheme for orthogonalization 
automatically switches from equation (38) to equation (43) across the concave 
region at the wing-body function and the values of the coefficient P are 
determined according to the near horizontality or vertically of the inner 
boundary. The resulting trajectories are smooth and swerve away from the 
orthogonal direction as they approach the outer boundary. A close-up view of 
the inner region is presented in Figure 32. Although the previous two examples 
involve analytically defined body surfaces, the procedure for generating 
orthogonal trajectories is equally applicable to body geometries input as 
discrete point sets, as is shown in the next example in Figure 33. This 
configuration contains sharp wing-body junctions and sharp wing tips, and 
orthogonality has been imposed only on the outboard portion of the wing. The 
resulting grid lines are smooth and no special problems are caused by the 
sharp wing tip as seen in the enlarged view in Figure 34.The next example in 
Figure 35 is a C-type grid for a section of a supersonic wing. The wing surface 
was input via discrete point coordinates. Orthogonality of grid lines near the 
leading edge is clearly visible in the enlarged view in Figure 36. 

As mentioned before, the use of the exponents P as defined by equations(42) 
and (47) ensures orthogonality, but not the prevention of grid lines of the same 
family from crossing each other. In body geometries involving sharp corners 
and concave regions the imposition of orthogonality at the surface causes grid 
lines to intersect within the grid field as shown in Figure 37. This problem occurs 
often in body geometries of realistic complexity and is usually untreatable by 
algebraic means. The homotopic procedure described here provides a way to 
maintain orthogonality at the inner boundary while preventing grid lines from 
crossing each other through the use of equation (48), where the coefficient q is 
defined by equation (49). The application of this procedure is demonstrated in 
the next few examples. The first example involves a blended wing-body 
configuration. A sectional grid generated with the orthogonality constraint 
imposed without the coefficient q results in the distorted grid shown in Figure 
38. The orthogonal trajectories emanating from the concave wing-body junction 
region are seen to have intersected those from the wing surface. The 
trajectories are heavily constrained to follow the nearly orthogonal direction 
until the last point immediately before the outer boundary causing distortion and 
depletion over a large portion of the grid. The use of equation (48) with q 
defined by 


q = 1.0 
k< 1.0 


N - 2 


-ik 



( 68 ) 


resolves the problem of intersecting trajectories and generates the smooth grid 
shown in Figure 39, with orthogonality preserved near the inner boundary. The 


coefficient K determines the rate of decay of orthogonality from the inner 
boundary to the outer boundary along trajectories. Values of K may be chosen 
to provide desired degrees of decay. The same principles are applied to the 
next case involving a blended wing-body-tail configuration. The basic 
orthogonlization scheme without the coefficient q, produces results similar to the 
previous example as shown in Figure 40. The strict adherence of the 
trajectories to the orthogonality constraint results in grid intersections and 
severe distortions over most of the surface. A smooth orthogonal grid shown in 
Figure 41 is obtained using equation (48) and (68) as in the previous example. 
It is important to note that any additional constraints, e.g., concentration of grid 
lines near the surface or grid clustering near a point on the surface are not 
altered by this orthogonalization procedure. 

All the previous grid examples were for external flow configurations where the 
grid was generated between the body surface as the inner boundary and a far- 
field surface as the outer boundary. The grid-generation procedure is equally 
applicable to some internal flow configurations e.g., ducts and nozzles as will 
be demonstrated in the following set of grid examples. In internal flow 
geometries all or part of the nozzle wall is usually chosen as the outer boundary 
while a number of choices exist for treating the inner boundary. These choices 
lead to three basic topological arrangements for internal grids. All three types 
are illustrated in the examples to follow. 

The first type of arrangement is very similar to a polar grid where the center-line 
of the nozzle is treated as the inner boundary and the nozzle wall serves as the 
outer boundary. A sectional grid with this arrangement is presented in Figure 
42. The inner boundary collapses to a point on the center-line in the planar 
section but the grid is generated by a procedure identical with that used for 
external flow configurations. Stacking of such planar grids results in a singular 
line along the center of the nozzle which is easily treatable in finite-volume 
applications. It may be noted that concentration and clustering of grids are still 
achievable by the methods used in external grids. 

The second type of grid arrangement has been applied to grid the exhaust from 
a blended wing-body aircraft in the region aft of the base as seen in Figure 43. 
An enlarged view is presented in Figure 44. In this arrangement the lateral 
diameter of the exhaust area is treated as the inner boundary and the upper 
and lower walls of the exhaust port is treated as outer boundaries. The grid is 
generated for the upper and lower halves of the region separately but complete 
continuity of the trajectories across the diameter is maintained. The trajectories 
are normal to the lateral diameter and the lines of the other family converge at 
the corner points fo the region at the ends of the diameter. In each sectional 
plane the trajectories of the internal grid match the external grid across the 
boundary of the internal grid. This is ensured by using the same distribution of 
points on this boundary when generating both the internal and external grids. A 
number of these sectional grids can be stacked behind the base of the body in 
order to study the flow in the wake region of the configuration. The blended 
wing-body configuration with four grids stacked in the wake region appears in 
Figure 45. 
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The third type of internal grids is different from the first two in that only a part of 
the duct wall is treated as the outer boundary. A sectional grid for a duct with 
circular cross-section is presented in Figure 46, in the left half plane. The grid is 
generated separately for the parts lying above and below the lateral diameter of 
the duct. The duct wall is considered in two parts. The part closer to the vertical 
diameter is chosen to be the outer boundary. The rest of the wall is treated as a 
side boundary. The distribution of points on the side boundary determines the 

distribution of the homotopy parameter n, used in blending. The Lateral diameter 
is treated as the inner boundary and continuity of the trajectories across it is 
maintained. This topological arrangement was employed to grid a duct with 
varying cross-sectional shapes. The duct shown in Figure 47, has cross- 
sections smoothly varying from a circle at one end to a narrow ellipse at the 
other end. The configuration with the sectional grids is presented in Figure 48. A 
smooth variation of the cross-sectional shapes ensures smooth grid lines in the 
longitudinal direction of the configuration and the resulting grid is suitable for 
computing internal flows in ducts. 

Application to solution of the Euler equations. 

In this section a sample set of results obtained by solving the Euler equations 
governing fluid flow is presented for a few configurations which were gridded by 
the procedure described in this report. All these results were obtained in the 
computational space discretized into hexahedral cells by a quasi-three- 
dimensional grid system created by stacking planar grids at stations along the 
longitudinal axis of the body. The Euler equations were solved on this 
discretized physical domain by a finite-volume algorithm with Runge-Kutta 
integration in time [1 1 ,12]. 

A plot of computed pressure contours in a sectional plane for a body with elliptic 

cross-sections at a Mach no of 1.6 and an angle of attack of 6° degrees is 
presented in Figure 49. The resulting contours are smooth and important 
details, e.g., the shock lying below the body are seen to be effectively captured. 
Surface pressure values for a cross section of an elliptic body are plotted 
against angular location in Figure 50. The Mach number and the angle of attack 
for this case are 2.5 and 5 degrees respectively. A reasonably close agreement 
is noted between the computed values and the experimental results given by 
Allen et.al [29]. A much more complex structure of shocks and expansions is 
seen to have been captured in the sectional plot of computed pressure contours 
for a blended wing body configuration presented in Figure 51. The flow 

conditions for this case were M = 3.0 and a = 3.0 degrees. Grids generated by 
the method described here are, therefore, seen to be suitable for use in 
computing high-speed flows for complex configurations. Computed pressure 
contours in the planform of a blended wing-body configuration appears in 
Figure 52. The shocks at the nose and the point where the wing begins to flare 
out are clearly seen to have been captured. The pressure contour lines are 
smooth in the longitudinal directions signifying that there is no severe distortion 
of grid cells in that direction. Finally, in Figure 53, a sectional plot of pressure 
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contours for a fighter aircraft-like configuration at M = 1.8 and a = 5° reveals 
all the expected details near the wing tip and the wing-body junction. 


Concluding remarks. 

Algebraic grid generation schemes are the fastest and the most efficient 
procedures available for generating computational grids for general 
configurations, because of the fact that the interpolant functions are specified 
explicitly. Acute sensitivity to boundary point distribution and the propagation of 
boundary discontinuities to points infield are, however, problems that have 
plagued algebraic schemes. Previously these difficulties were overcome by 
clustering the boundary mesh points at points with slope discontinuity and local 
controls on the interpolants e.g., in multisurface transformations. These local 
controls are very much dependent on the particular configuration under 
consideration and apriori knowledge of the problem areas are usually 
necessary in formulating the control functions. In the present homotopic 
procedure the interpolants are modified by control of the homotopy parameter 
where the control functions are derived from the boundary data and play roles 
analogous to control functions in elliptic systems. This homotopic method 
prevents propagation of discontinuities and intersection of the trajectories while 
maintaining near orthogonality at the boundary. The resulting grids are smooth 
and have been successfully used in high-speed Euler computations. 

All grid generation schemes, however complex, need some degree of human 
intervention in order to achieve grids with desired characteristics for sufficiently 
complex configurations. The human intervention is best implemented in 
interactive computer graphics environments; and algebraic schemes, by virtue 
of their speed, are ideally suited for this kind of applications. For a quasi-three- 
dimensional grid system, each planar grid can be viewed on a terminal and 
modified by changing the grid parameters if desired in a short amount of time 
provided the communication rate is sufficiently high. Problem areas can be 
identified by viewing grids with a small number of points thus further reducing 
development time. 
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Figure 2 


Schematic diagram of quasi-three-dimensional grid 
systems 
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Figure 5.b Family of curves after perturbation. 

Figure 5 Homotopic perturbation to achieve near-orthogonality. 





Figure 6.b Region of intersection. 



Figure 6.c Corresponding intersecting maps. 


Figure 6 Intersection of orthogonal trajectories. 
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Figure 7.a 

Modification of the x-coordinate for horizontal boundaries. 



Figure 7.b Modification of the y-coordinate for vertical boundaries. 

Figure 7 Modification of coordinates to achieve orthogonality at 

horizontal and vertical boundaries. 
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Figure 8 


Vectors used in orthogonalization process. 
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Figure 9 


Body geometry defined as transition surface. 
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Figure lO.a 

Normalized nose section at X = Xpg. 




Figure lO.b 

Normalized terminal section at X = X 


Figure 10 Normalized end shapes. 



Figure 11 


Normalized shape of grid outer boundary. 
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Figure 21 .b Cut-away view of three-dimensional grid. 
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Figure 21 .c Fighter aircraft with two planar grids. 
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Figure 25 Grid lines clustered at tip of elliptic body. 
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Figure 42 Polar-like grid for nozzle with grid clustering. 
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Figure 43 Grids in exhaust area behind base of blended wing-body 

configuration. 




Figure 44 Continuity between inner and outer grids in exhaust area. 
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Figure 48 Duct configuration with sectional grids. 
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Figure 49 Computed pressure contours for elliptic body M = 1.6, 

a = 6°. 


Surface pressure on elliptic body 

9 = 0 ° 



Computed surface pressure on elliptic body. 0° ^0 ^ 180 
M = 2.5, a = 5°. 


Figure 50 






